{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sisl\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example we will begin to cover some of the extraction utilities in `sisl` allowing one to really go *in-depth* on analysis of calculations.\n",
    "\n",
    "We will begin by creating a large graphene flake, then subsequently a hole will be created by removing a circular shape.\n",
    "\n",
    "After the calculations, you are encouraged to explore the `sisl` toolbox for ways to extract important information regarding your system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "graphene = sisl.geom.graphene(orthogonal=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "elec = graphene.tile(25, axis=0)\n",
    "H = sisl.Hamiltonian(elec)\n",
    "H.construct(([0.1, 1.43], [0., -2.7]))\n",
    "H.write('ELEC.nc')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "device = elec.tile(15, axis=1)\n",
    "center = device.center(what=\"xyz\")\n",
    "device = device.remove(\n",
    "    device.close(center, R=10.)\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will also (for physical reasons) remove all dangling bonds, and secondly we will create a list of atomic indices which corresponds to the atoms at the edge of the hole. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dangling = [ia for ia in device.close(center, R=14.)\n",
    "                if len(device.close(ia, R=1.43)) < 3]\n",
    "device = device.remove(dangling)\n",
    "edge = []\n",
    "for ia in device.close(center, R=14.):\n",
    "    if len(device.close(ia, R=1.43)) < 4:\n",
    "        edge.append(ia)\n",
    "edge = np.array(edge)\n",
    "\n",
    "# Pretty-print the list of atoms (for use in sdata)\n",
    "# Note we add 1 to get fortran indices\n",
    "print(sisl.utils.list2str(edge + 1))\n",
    "\n",
    "Hdev = sisl.Hamiltonian(device)\n",
    "Hdev.construct(([0.1, 1.43], [0, -2.7]))\n",
    "Hdev.geometry.write('device.xyz')\n",
    "Hdev.write('DEVICE.nc')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercises\n",
    "\n",
    "Please carefully go through the `RUN.fdf` file before running TBtrans, determine what each flag means and what it tells TBtrans to output.  \n",
    "Now run TBtrans:\n",
    "\n",
    "    tbtrans RUN.fdf\n",
    "    \n",
    "In addition to the previous example, many more files are now being created (for all files, the `siesta.TBT.AV*` files are the $k$-averaged equivalents). You should, while reading the below list, also be able to specify which of the fdf-flags are responsible for the creation of which files.\n",
    "\n",
    "- `siesta.TBT.ADOS_*`  \n",
    "   The $k$-resolved spectral density of states injected from the named electrode\n",
    "- `siesta.TBT.BDOS_*`  \n",
    "   The $k$-resolved bulk Green function density of states for the named electrode\n",
    "- `siesta.TBT.BTRANS_*`  \n",
    "   The $k$-resolved bulk transmission through the named electrode\n",
    "- `siesta.TBT.DOS`  \n",
    "   Green function density of states\n",
    "- `siesta.TBT.TEIG_<1>-<2>`  \n",
    "   Transmission eigenvalues for electrons emitted from `<1>` and injected in `<2>`.\n",
    "   \n",
    "\n",
    "This exercise mainly shows a variety of methods useful for extracting data from the `*.TBT.nc` files in a simple and consistent manner. You are encouraged to play with the routines, because the next example forces you to utilize them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tbt = sisl.get_sile('siesta.TBT.nc')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As this system is not a pristine periodic system, we have a variety of options available for analysis.\n",
    "\n",
    "First and foremost, we plot the transmission:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(tbt.E, tbt.transmission(),label='Av');\n",
    "plt.plot(tbt.E, tbt.transmission(kavg=tbt.kindex([0,0,0])), label=r'$\\Gamma$'); plt.legend()\n",
    "plt.ylabel('Transmission'); plt.xlabel('Energy [eV]'); plt.ylim([0, None]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "The contained data in the `*.TBT.nc` file is very much dependent on the flags in the fdf-file. To ease the overview of the available output and what is contained in the file, one can execute the following block to see the content of the file.\n",
    "Check whether the *bulk transmission* is present in the output file and if so, add it to the plot above to compare the bulk transmission with the device transmission.  \n",
    "\n",
    "There are two electrodes, hence two bulk transmissions. Is there a difference between the two? If so, why, if not, why not?\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(tbt.info())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Density of states"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "source": [
    "We may also plot the Green function density of states as well as the spectral density of states from each electrode:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DOS_all = [tbt.DOS(), tbt.ADOS(0), tbt.ADOS(1)]\n",
    "plt.plot(tbt.E, DOS_all[0], label='G');\n",
    "plt.plot(tbt.E, DOS_all[1], label=r'$A_L$');\n",
    "plt.plot(tbt.E, DOS_all[2], label=r'$A_R$');\n",
    "plt.ylim([0, None]); plt.ylabel('Total DOS [1/eV]'); plt.xlabel('Energy [eV]'); plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Can you from the above three quantities determine whether there are any localized states in the examined system?  \n",
    "*HINT*: What is the sum of the spectral density of states ($\\mathbf A_i$) compared to the Green function ($\\mathbf G$) density of states?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Examining DOS on individual atoms\n",
    "\n",
    "The total DOS is a measure of the DOS spread out in the entire device region. However, TBtrans calculates, and stores all DOS related quantities in orbital resolution. I.e. we are able to post-process the DOS and examine the atom and/or orbital resolved DOS.  \n",
    "To do this the `.DOS` and `.ADOS` routines have two important keywords, 1) `atom` and 2) `orbital` which may be used to extract a subset of the DOS quantities. For details on extracting these subset quantities, please read the documentation by executing the following line:\n",
    "\n",
    "    help(tbt.DOS)\n",
    "    \n",
    "The following code will extract the DOS _only_ on the atoms in the hole edge."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DOS_edge = [tbt.DOS(atoms=edge), tbt.ADOS(0, atoms=edge), tbt.ADOS(1, atoms=edge)]\n",
    "plt.plot(tbt.E, DOS_edge[0], label='G');\n",
    "plt.plot(tbt.E, DOS_edge[1], label=r'$A_L$');\n",
    "plt.plot(tbt.E, DOS_edge[2], label=r'$A_R$');\n",
    "plt.ylim([0, None]); plt.ylabel('DOS on edge atoms [1/eV]'); plt.xlabel('Energy [eV]'); plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Comparing the two previous figures leaves us with little knowledge of the DOS ratio. I.e. both plots show the _total_ DOS, and they are summed over a different number of atoms (or orbitals if you wish). Instead of showing the total DOS, we can normalize the DOS by the number of atoms; $1/N_a$ where $N_a$ is the number of atoms in the selected DOS region. With this normalization, we can compare the average DOS on all atoms with the average DOS on only edge atoms.  \n",
    "\n",
    "---\n",
    "\n",
    "The `tbtncSileTBtrans` can readily do this for you.  \n",
    "Read about the `norm` keyword in `.DOS`, and also look at the documentation for the `.norm` function:\n",
    "\n",
    "    help(tbt.DOS)\n",
    "    help(tbt.norm)\n",
    "\n",
    "---\n",
    "\n",
    "Now create a plot with DOS normalized per atom by editing the below lines, feel free to add the remaining DOS plots to have them all:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_all = tbt.norm(<change here>)\n",
    "N_edge = tbt.norm(<change here>)\n",
    "plt.plot(tbt.E, DOS_all[0] / N_all, label=r'$G$');\n",
    "plt.plot(tbt.E, DOS_edge[0] / N_edge, label=r'$G_E$');\n",
    "plt.ylim([0, None]); plt.ylabel('DOS [1/eV/atom]'); plt.xlabel('Energy [eV]'); plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### DOS depending on the distance from the hole\n",
    "\n",
    "We can further analyse the DOS evolution for atoms farther and farther from the hole.  \n",
    "The following DOS analysis will extract DOS (from $\\mathbf G$) for the edge atoms, then for the nearest neighbours to the edge atoms, and for the next-nearest neighbours to the edge atoms.  \n",
    "Try and extend the plot to contain the DOS of the next-next-nearest neighbours to the edge atoms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get nearest neighbours to the edge atoms\n",
    "n_edge = Hdev.edges(edge, exclude=edge)\n",
    "# Get next-nearest neighbours to the edge atoms\n",
    "nn_edge = Hdev.edges(n_edge, exclude=np.hstack((edge, n_edge)))\n",
    "# Try and create the next-next-nearest neighbours to the edge atoms and add it to the plot\n",
    "plt.plot(tbt.E, tbt.DOS(atoms=edge, norm='atom'), label='edge: G');\n",
    "plt.plot(tbt.E, tbt.DOS(atoms=n_edge, norm='atom'), label='n-edge: G');\n",
    "plt.plot(tbt.E, tbt.DOS(atoms=nn_edge, norm='atom'), label='nn-edge: G');\n",
    "plt.ylim([0, None]); plt.ylabel('DOS [1/eV/atom]'); plt.xlabel('Energy [eV]'); plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Learned methods\n",
    "\n",
    "- fdf-flags for TBtrans to specify which quantities to calculate\n",
    "- Opening files via `get_sile`\n",
    "- Extraction of DOS with different normalizations.\n",
    "- Extraction of DOS for a subset of atoms/orbitals.\n",
    "- Determining coupling atoms "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
